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It is experimentally observed and theoretically proved that the distance between topmost layers of 
a metal surface has a contraction. However, well-known potentials such as Lennard- Jones and Morse 
potentials lead to an expansion of the surface inter-layer distance. Such simple potentials therefore 
cannot be used to study metal surface relaxation. In this paper, extensive Monte Carlo simulations 
are used to study the silver (111) surface with both the Gupta potential (GP) and the Embedded 
Atom Method (EAM) potential. Our results of the lattice relaxation at the (111) surface of silver 
show indeed a contraction for both potentials at low temperatures in agreement with experiments 
and early theories. However at higher temperatures, the EAM potential yields a surface melting at 
~ 700 K very low with respect to the experimental bulk melting at ~ 1235 K while the GP yields a 
surface melting at ~ 1000 K closer to the bulk one. In addition, we observe with the EAM potential 
an anomalous thermal expansion, i. e. the surface contraction becomes a surface dilatation with 
respect to the bulk, at ~ 900 K. The Gupta potential does not show this behavior. We compare our 
results with different experimental and numerical results. 
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I. INTRODUCTION 

It is well known that surface atoms of a material have 
a behavior different from that of bulk atoms, mainly be- 
cause of the lack of neighbors and the surface geometry. 
Different kinds of surface behavior can occur according 
to the nature of the material and the surface orientation. 
The case of metallic materials has been well studied the- 
oretically [1^ 71 and by means of different experimental 
techniques [8l4lOl|. 

The contraction of the lattice spacing is experimentally 
observed by Medium Energy Ion Backscattering (MEIS) 
3, Low- Energy Electron Diffraction (LEED) [9] or X-ray 
scattering . Other methods such as elastic He scatter- 
ing and electron energy-loss spectroscopy have also been 
used. In a theoretical point of view, Gupta has shown an- 
alytically that for classical pairwise potentials (Lennard- 
Jones, Morse, ...), inter-layer distance near the surface 
have an expansion. He has also shown that the Tight- 
Binding-Potential (TBP) and the so-called Gupta poten- 
tial (GP) lead to a contraction of inter-layer distance at 
metallic surfaces. 

Despite the important number of experimental and 
theoretical techniques used to observe this phenomenon, 
there exists an unsolved question on how the contrac- 
tion evolves with increasing temperature. There are two 
contradictory answers in the literature: X-ray scatter- 
ing 11] and LEED "9] as well as molecular dynamics 
(MD) simulations using the Embedded Atom Method 
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(EAM) [3, show the surface inter-layer distance al- 
ways smaller than the bulk one at the same temperature, 
namely surface is contracted, whereas MEIS experiments 
[§] and ab-initio density-functional theory (DFT) calcu- 
lations [13] show an anomalous thermal expansion of the 
surface at some temperature below the bulk melt. Facing 
this long-standing unsolved question, we wanted to carry 
out a Monte Carlo (MC) study in an attempt to answer 
that question. To our knowledge, there is no MC simu- 
lations in literature so far about this subject, although 
some MC simulations have been used to reproduce exper- 
imental patterns such as surface blocking pattern during 
scattering process [see Ref. [13] for the Pb (110) surface]. 
Given a tremendous number of numerical studies on sur- 
face problems, it is surprising that no MC simulation has 
been performed so far to see the variation of the surface 
intcr-layer distance. 

The purpose of this paper is thus to investigate by MC 
simulation the variation of the lattice spacing between 
the topmost layers of the (111) silver surface versus tem- 
perature. 

In order to simulate such a behavior as accurately as 
possible, we have considered potentials which describe as 
well as possible the material. The EAM potential is often 
used in MD simulations and especially for the Ag (111) 
surface 0, [13: [31 ■ Working with this potential allows us 
to compare our results with other numerical studies using 
the same potential. Furthermore, the EAM potential re- 
produces accurately the bulk melting temperature of Ag. 
On the other hand, the Gupta potential describes well 
surface and cluster behaviors [l5l - [l9| . The melting tem- 
perature of bulk Ag is also well reproduced with this po- 
tential. That was the reason why the two potentials GP 
and EAM have been used since many years to simulate 
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silver material and other metallic crystals. However, as 
will be seen below, the two potentials, although yielding 
the same result for low-temperature surface contraction, 
give different results at higher temperature concerning 
the surface contraction and the surface melting. 

In Section [n] we recall essential properties of the two 
potentials with their sets of parameters and we briefly de- 
scribe our algorithm. The results obtained by our com- 
putation are shown in Section IIIIl Concluding remarks 
are given in Section |IVl 



II. MODEL AND MONTE CARLO METHOD 
A. Potentials 

1. Gupta potential 

In computer simulation of metals, the Gupta potential 
is one of the most used semi-empirical potentials. There 
is a multiple reason for this success but the main ones 
are the accuracy of its results for metals and its easy and 
quick implementation. The Gupta potential is derived 
from the Gupta's expression of the cohesive energy of 
the bulk material and is based on the second-moment 
approximation of the electron density of states in the 
tight-binding theory. It includes implicitly some many- 
body interactions. The expression of the potential is the 
following: 
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where A is a constant given in eV, Tq the equilibrium 
nearest-neighbor (NN) distance in the bulk metal, p the 
repulsive interaction range and q the attractive one. j 
is the distance between the atoms i and j and £"„ an 
energy constant which depends on the size of the system. 
In Eq. ([T|), the first part gives is a Born-Mayer pairwise 
repulsion energy term and the second part is the many- 
body attractive contribution. 

The parameters used in this work are reported in Table 

III 



Table I. Gupta parameters for silver 



Parameter 


Value'' 


A (eV) 


0.09944 


P 


10.12 


g 


3.37 


ro(A) 


2.892 


En 


2.52_^ 



" Reference Il6l 
Fitted with experimental bulk melting temperature in [20 



2. EAM potential 



Several authors |21| - |23l | have proposed a method based 
on density-functional ideas called EAM. We used a ver- 
sion of the EAM potential given in Ref. [2^ . The param- 
eters of the EAM potential used in this work for silver 
are reported in Table HIl In the EAM potential, the total 
potential energy is given by: 



(2) 



where (j^ij represents the pair energy between two atoms 
i and j at the distance r^j . Fi is the embedding energy 
function which represents the energy to embed an atom 
into a local site with electron density pi. The electron 
density pi has the following expression: 
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where r^., /e, /3 and A are constant parameters which 
are given in Table HH The embedding function has the 
following form: 



Fip)={ EtoF^ " ^■^^P'^ 
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The pair energy expression between atoms i and j is : 
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All constant parameters in the above expressions such as 
F„. etc. are listed in Table HIl Note that MD simulations 
using EAM potential yield the bulk melting temperature 
at 1170 K for Ag [22| which is to be compared to the 
experimental value of 1235 K [23 |. 



B. The algorithm 

The c omp lete description of the algorithm can be found 
in Ref. [2^. For the present surface problem, we have 
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Table II. EAM parameters for silver 



Parameter 


Value'" 


Parameter 


Value 


• e 


2.891814 




-0.221025 








541558 


pe 


15.539255 


p 

77.3 


-0.967036 


Q 


7.944536 


Fo 


-1.75 




4.237086 


Fi 





A 


0.266074 


F2 


0.983967 


B 


0.386272 


F3 


0.520904 


K 


0.425351 


V 


1.149461 


A 


0.850703 




-1.751274 




-1.729619 







" Reference |23| 



used the the periodic boundary conditions (PBC) in two 
directions x and y. Let us briefly recall the main steps of 
the algorithm. The step that requires a long CPU time is 
the construction of the PBC table which contains neigh- 
bors of each atom of the simulation sample. The reader 
is referred to Ref. [1^ for details. Note that to apply 
correctly the PBC we have to take into account instan- 
taneous fluctuations of the sample size in each direction 
at every MC step. 

Using the least square approximation for experimen- 
tal data of the NN distance measured for a bulk sample 
taken from Ref. [l^l we obtain the linear thermal ex- 
pansion (i7VAr(A) = ulT + 2.87053 with the coefficient 
aL = 6.27617 x 10^^ . In Fig. [T]we have plotted exper- 
imental data taken from Ref. [24| and the least square 
approximation used in this paper. The lattice constant 
at a given temperature is now an input data. In our 
calculations, we begin at T = 200 K in order to stay in 
the classical limit. Indeed, as we can see in Ref. [2^ 
that at lower T, the thermal expansion coefficient does 
not have the same linear behavior below and above 200 
K. At each temperature, we start the simulation with 
the lattice constant taken from Fig. [T] In doing so, we 
favor a fast convergence to equilibrium. This was what 
we observed in test runs. Only the surface atoms and 
an adjustable number of layers beneath the surface are 
moving. It allows us to choose a number of moving lay- 
ers to produce precise results while saving CPU time. Of 
course, we have to make sure that final results do not de- 
pend on the choice of the number of moving atoms (see 
an example in Fig. [5]). 

Note that by using the correct value of the lattice con- 
stant at each temperature for interior layers far from the 
surface, we take into account some temperature effects 
but we neglect local fiuctuations in those far layers. The 
effect of such approximation has been tested: the result 
confirms that this does not affect our conclusion. As 
seen below, in most simulations which have been made, 
atoms are allowed to move in the first three layers. We 
have observed that the contraction of the surface lattice 
spacing disappears rapidly at a few layers from the sur- 
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Figure 1. (Color online) Experimental data by X-rays diffrac- 
tion measurements(green circles) from Ref. [2J] and the least 
square linear fit (red line). 



face. This approximation is justified by tests we have 
made and by experimental observations and theoretical 
calculations mentioned above. 




Figure 2. (Color online) Side view of a silver sample, with 
periodic boundary conditions, after several millions of MC 
steps at 300 K. Only atoms of the three topmost layers are 
allowed to move. See text for explanation. 

Note that we use the Metropolis updating criterion. 
At a given temperature, we equilibrate the system before 
computing the averages of quantities needed for analyz- 
ing the system behavior as functions of temperature. 



III. RESULTS AND DISCUSSIONS 

In this section we describe our MC results for the (111) 
surface of silver. We have studied a system of 256 atoms 
per layer in which atoms in the first three layers are al- 
lowed to relax at each MC step. From the fourth layer 
inward, the atoms dilate uniformly with temperature us- 
ing data of Fig. [T] For comparison, we have also con- 
sidered a system of 100 atoms per layer with 3 moving 
layers. As it turned out, the thermal relaxation of the 
first moving layers is not significantly different for these 
two surface sizes. Therefore, we show below our results 
obtained with the surface size of 256 atoms on the top 
of five beneath layers (6 layers in all) with the PBC in 
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the X and y directions. Usually, our simulations are car- 
ried out over 5 to 6 millions of MC steps per atoms. We 
discarded the two first millions during equilibrating time 
and averaged physical quantities over the following 2-3 
millions of MC steps. 



A. Surface melting 

The first step is the monitoring of the surface stabil- 
ity versus temperature. Of course, our calculations have 
been carried out for temperatures below the bulk melt- 
ing temperature (Tm) because we use the temperature- 
dependent bulk lattice spacing below Tm shown in Fig. [1] 
for bulk atoms. Note however that very close to r„i, the 
thermal coefficient does not follow the linear behavior. 
We have therefore heated the system from the ground 
state to temperatures between 200 K and 1200 K, just 
below Tjn, with a step of 50 K. 

Distance (A) 




1000 1200 1400 
T (K) 



Figure 3. (Color online) Mean displacement amplitude (in A) 
between atoms and their perfect ground-state lattice positions 
computed with Gupta potential. 

The mean displacement amplitude allows us to estab- 
lish the validity domain of temperatures. This quantity 
is computed as follows : at each MC step, we compute a 
spatial average of all separation distances between atoms 
with their ground-state node positions and then we com- 
pute its thermal average through MC steps. As we can 
see in Fig. [31 with the Gupta potential, after 1000 K the 
atoms undergo a sudden jump indicating a phase change: 
the surface becomes disordered (melted) between 1000 
K and 1050 K. This value from our MC simulation us- 
ing Gupta potential is very close to the value 1100 K 
obtained from MD simulation using the EAM potential 
[12]. A study is needed to determine the precise surface 
melting temperature but it is not the aim to search for 
such a precision of this paper. 

For the EAM potential, we have obtained a result dif- 
ferent from that of MD simulation: the surface melting 
temperature takes place at about 700 K, while MD sim- 
ulation fonds it at 1100 K as we mentioned above. In 
order to see clearly this surface melting, we have com- 



puted with the EAM potential the structure factor for 
the topmost layers. The structure factor is computed 
as follow: 



Ni 



(7) 



where dj is the position vector of an atom in the layer, Ni 

the number of atoms in a layer and the reciprocal lat- 
tice vector which has the following coordinate (in reduced 
units): 27r (l; — VS; O) . The angular brackets < ... > in- 
dicate thermal average taken over MC run time. The 
above " order parameter" , which allows us to monitor the 
long-range surface order, is plotted for the surface layer 
in Fig. m As we can see, the long-range order is lost 
at ~ 700 K. In order to investigate in more details the 
surface melting, we have also computed the Oq order pa- 
rameter introduced in Ref. [l^]. This order parameter 
describes the short-range hexatic orientational order of 
the surface which is defined as follows: 
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with 
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(9) 



where the sum runs over the NN pairs and Qjk is the an- 
gle which the j — k bond, when projected on the xy plane, 
forms with the x axis. The 5 parameter is taken as one- 
half the average inter-layer spacing. The weighting func- 
tion, Wjk, allows us to differentiate the "non coplanar" 
and the "coplanar" neighbors. With a coplanar neigh- 
bor, the weighting function takes a maximum value. We 
have to calculate the spatial average of Oq taken over all 
atoms of the surface layer and then calculate its thermal 
average over MC run time. We plot the so-averaged Oq 
parameter and the Sj^ structure factor versus tempera- 
ture in Fig. m 

We observe that the long-range and short-range orders 
are lost at the same temperature: the melting of the sur- 
face occurs at around 700 K where both order parameter 
change their curvature. 

Snapshots of the system at two temperatures close to 
the surface melting are shown in Figs. [5] and El for the 
two potentials. 

The difference in the surface melting temperature for 
the two potentials is surprising since both give the same 
bulk melting temperature. One possible explanation 
comes from the construction of their set of parameters: 
the Gupta potential was fitted in Ref. [2^ with the bulk 
melting temperature, while the EAM potential is fitted to 
reproduce physical quantities such as the lattice constant, 
the elastic constants Cii,Ci2 and C44 and sublimation 
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Figure 4. (Color online) Structure factor Sjf (green triangles) 
and Oe order parameter (red circles) of the first layer versus 
temperature for the EAM potential. 



B. Surface contraction 

In order to see the variation of the inter-layer distance 
at the surface, we have plotted in Fig. [7] the z-position 
distribution obtained for the two potentials EAM and 
GP. As we can see, when the temperature increases, the 
three peaks which represent the z coordinates of the three 
moving layers, are shifted to the left, indicating a con- 
traction toward the bulk. Note that the two potentials 
lead to the same behavior. 

In these two z-axis distribution functions, we can see 
that at a given T, the deviation of the topmost layer from 
its perfect position is largest. This deviation becomes 
smaller for the next two layers and it should disappear 
at a few layers from the surface. This observation justifies 
our approximation to let only the three topmost layers 
to relax. Previous theoretical results and experimental 
data also give support to our hypothesis. 



energy, etc. but not with the bulk melting temperature. 
A consequence is that surface atoms in the Gupta case 
are more strongly attached to the interior part, so that 
we need a much higher temperature to make the sur- 
face melt. We believe that the high "stability" of surface 
atoms is the reason why the Gupta potential was very 
suitable for the study of clusters of very few number of 
atoms [Tsj . 

For both potentials, a contraction of the lattice spacing 
between two topmost layers occurs but its temperature 
dependence is different in the two cases, as will be shown 
below. 



z-distribution 



EAM 





Figure 5. Pictures of the system for two different tempera- 
tures computed witli the EAM potential: 550 K (left) and 
650 K (right). Only three layers are allowed to relax. 
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Figure 6. Pictures of the system for two different tempera- 
tures computed with the Gupta potential : 950 K (left) and 
1050 K (right). Only three layers are allowed to relax;. 



Figure 7. (Color online) z-position distribution at different 
temperatures calculated with EAM and Gupta potentials. 
The surface at its non relaxed position corresponds to z—Q. 

In order to examine closely the inter-layer contraction, 
we have plotted in Fig. [8] the inter-layer distance A12 
between layer 1 and 2, at different temperatures. This 
quantity is determined by averaging over the correspond- 
ing peaks during MC simulation time. We only focus our 
attention on A12 since for this quantity there is a gen- 
eral agreement between experiments and between exper- 
iments and theories. For inner layers, there is no such 
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agreement between experiments and theories, and even 
between experiments, depending on the technique used: 
LEED 0, MEIS [3, ... 

Let us comment the results shown in Fig. |8l Several 
remarks are in order: 

(i) At low temperatures, for both potentials, A12 has 
a contraction of about 2.5% at 300 K, compared to a 
contraction of 0.5% with MD simulations using EAM 
potential (JH^l- Experimental contraction of 2.5% was 
observed in Ref. Q. 

(ii) The distance A12 increases with increasing tem- 
perature: in the case of EAM potential, A12 crosses the 
bulk limit at ~ 900 K, indicating a surface expansion, 
in agreement with experimental results in Ref. [8]. This 
"expansion" was not found by MD simulations with EAM 

otential (see Refs. 0,[ill) up to 1100 K and in LEED 
A surface expansion was observed in an ab-initio 
DFT calculation 

(iii) As we can see in Fig. El MC results with GP does 
not show the above-mentioned anomalous expansion up 
to 1150 K slightly above the surface melting temperature 
observed above with this potential. We will give below 
an explanation about this point at least from our MC 
results. 
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and MD result using the same potential resides not only 
in the difference of surface melting temperatures (ours is 
700 K, the MD one is 1100 K) but also in the existence 
of the anomalous expansion. 

For GP, we do not observe the anomalous expansion 
up to 11500 K as shown in Fig. |S1 So the existence or not 
of an anomalous expansion depends on the potential at 
least with our MC results. Let us explain this as follows: 

(i) With the EAM potential: The surface melts at a 
temperature much lower than that of the anomalous ex- 
pansion. As we can see in Fig. [51 at around 700 K, the 
melting of the first layer causes the apparition of a lit- 
tle peak but A12 is still smaller than the bulk inter-layer 
distance. The fact that the anomalous expansion does 
not occur at the melting temperature of the first layer 
means that the melted surface is in a "two-dimensional" 
liquid state in the temperature region between 700 K and 
900 K. This liquid surface layer is "detached" from the 
remaining crystal only at temperatures higher than 900 
K. More accurate experiments of surface melting are de- 
sirable to allow a better understanding of this point. 

(ii) With the GP: Since the surface melts at ~ 1000 — 
1050 K (see Fig. [3]), if there is an anomalous expansion, 
this should happen at a higher temperature. Given the 
fact that with EAM, our MC result indicates the anoma- 
lous expansion occurring 200 K above the surface melting 
temperature, we can imagine the same scenario with the 
GP: an expansion may be possible at ~ 1200 — 1235 K, 
namely at the bulk melting. 

To close this section, let us summarize here that (i) 
the surface melting occurs below the bulk one, but the 
distance to the bulk melting temperature depends on the 
potential: the EAM and Gupta potentials give the same 
bulk melting temperature but different surface melting 
temperatures (ii) an anomalous dilatation of surface is 
possible only at a temperature much higher than the sur- 
face melting temperature. 



Figure 8. (Color online) Inter-layer A12 computed with 
Gupta and EAM potentials. Dashed line represents the bulk 
inter-layer distance. See text for comments. 

As said above, the contraction of the first inter-layer 
distance A12 has been observed by both MEIS and 
LEED at low temperatures. For A23 , the distance be- 
tween the second and third layers, these two methods are 
not in agreement: at low T, MEIS finds a little expan- 
sion of A23 while LEED finds a contraction. Here, our 
computations are in agreement with LEED results. 

Let us discuss about the controversial point on the 
anomalous expansion of A12, namely the transition from 
contraction to expansion when A12 crosses the bulk line 
(see Fig. [5]). As said earlier, this has been experimentally 
observed by MEIS at about 750 K Q. We have found 
this with our simulations using EAM potential at ~ 900 
K. Our expansion is about 5% at 1150 K, with respect to 
the bulk value. The disagreement between our MC result 



IV. CONCLUDING REMARKS 

Monte Carlo simulations with both EAM and Gupta 
potentials show two different surface melting tempera- 
tures for the (111) surface of a silver sample: while both 
EAM and Gupta potentials give the same bulk melting 
temperature (~ 1200 K), the EAM potential yields a sur- 
face melting at ~ 700 K and the Gupta potential shows 
a surface melting at ~ 1000 — 1050 K close to the bulk 
melting. 

However, both potentials show a contraction of the 
topmost inter-layer distance of that surface at low tem- 
peratures, in good agreement with experiments and the- 
ories. Our results of the temperature-dependence of 
the inter-layer contraction indicate a strong potential- 
dependence: the variation of the contraction with in- 
creasing temperature shows a difference of the two po- 
tentials. Using the EAM potential, our MC simulations 
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show an anomalous thermal expansion, namely the dis- 
tance between the topmost layers is larger than that be- 
tween two adjacent bulk layers. This can be explained 
by the fact that the anomalous thermal expansion oc- 
curs only at a temperature much higher than the sur- 
face melting temperature: the already-disordered surface 
layer is loosely attached to the crystal. This surface "di- 
latation" is in agreement with MEIS experiments and 
ab-initio DFT calculations. On the other hand, using 
the Gupta potential, our results show that the surface 
inter-layer distance varies with the temperature but it is 
always smaller than the bulk distance at least up to tem- 
peratures close to the bulk melting. We attribute the non 
observation of a surface dilatation in this case by the fact 



that surface melting occurs too close to the bulk melting. 
The high surface melting temperature has been found in 
MD simulations with EAM potential [l^. The existence 
of surface dilatation has been observed in MEIS experi- 
ments |8i but not in LEED experiments [9] and in X-ray 
scattering We think that this difference comes in 

part from the difference of surface flatness, surface clean- 
ness (contaminated or not) and experimental conditions. 
We hope that experimentalists could resolve this impor- 
tant question on anomalous high-temperature behavior 
of Ag(lll). Experiments to determine the melting tem- 
perature of the Ag (111) surface is also needed to choose 
a potential appropriate for Ag crystal. 
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